#remotes::install_github("thomasp85/patchwork")
#remotes::install_github("tidyverse/ggplot2", ref = remotes::github_pull("5592"))
#library(optparse)
#library(readr) tidyverse
#library(stringr) tidyverse

# R 4.3.2

#if (!require("BiocManager", quietly = TRUE))
#  install.packages("BiocManager")

library(tidyverse)
library(GGally)
library(ggrepel)

library(DESeq2)
library(edgeR)      # cpm
library(vsn)        # meanSdplot
library(ReportingTools)  # HTMLReport
library(ape)
library(ShortRead)
library(apeglm)
library(EDASeq)  # plotRLE

library(mixOmics)
library(matrixStats)

library(RColorBrewer)
library(viridis)
library(patchwork)
library(ggpointdensity)
library(data.table)
library(DT)
library(htmlTable)

library(DGEobj.utils)

source("/users/jaritz/Gaidt/bioinf/software/scripts/deseq2.plot.R")
source("/users/jaritz/Gaidt/bioinf/software/scripts/multiple_matrix_DESeq2.functions.R")  # mPCA
source("/users/jaritz/Gaidt/bioinf/software/scripts/analyse.functions.R")  # write_extended_results_table
workdate <- "2025Oct13"
workdir <-  paste0("/groups/gaidt/user/markus.jaritz/projects/Markus/",workdate,"/")

sampleTable=paste0(workdir,"doc/samples.txt")

dir.create(paste0(workdir,"/results"),recursive = T,showWarnings = F)
setwd(paste0(workdir,"/results"))
getwd()

odir <- getwd()

params=list(
  sampleTable='/users/jaritz/Gaidt/bioinf/software/danbing-tk/samples_experiment.txt',
  conditions_colname='condition',
  control_condition='exp1',
  experiment_condition='exp2',
  countTableAnnot='/scratch-cbe/users/markus.jaritz/Gaidt/Moritz/2025Oct13/results/featureCount_concat/featureCount_counts.annot.txt',
  countTable='/scratch-cbe/users/markus.jaritz/Gaidt/Moritz/2025Oct13/results/featureCount_concat/featureCount_counts.txt',
  odir=odir,
  downsample_pairs=0,
  downsample_ma=0,
  session="http://localhost:60151/goto?"
)
params
p_ <- GGally::print_if_interactive
set.seed(42)
p_ <- GGally::print_if_interactive
set.seed(42)

print(paste("sampleTable: ",params$sampleTable))
[1] "sampleTable:  samples_experiment.txt"
print(paste("conditions_colname: ",params$conditions_colname))
[1] "conditions_colname:  condition"
print(paste("control_condition: ",params$control_condition))
[1] "control_condition:  exp1"
print(paste("experiment_condition: ",params$experiment_condition))
[1] "experiment_condition:  exp2"
print(paste("countTable: ",params$countTable))
[1] "countTable:  /groups/gaidt/user/markus.jaritz/projects/Markus/2025Oct13/results/featureCount_concat/featureCount_counts.txt"
print(paste("countTableAnnot: ",params$countTableAnnot))
[1] "countTableAnnot:  /groups/gaidt/user/markus.jaritz/projects/Markus/2025Oct13/results/featureCount_concat/featureCount_counts.annot.txt"
print(paste("downsample_pairs: ",params$downsample_pairs))
[1] "downsample_pairs:  1000"
print(paste("downsample_ma: ",params$downsample_ma))
[1] "downsample_ma:  0"
# DEBUG
# params <- list(celltype="BLaER1")

1 INPUT

#
# count table
#
t <- read_delim(params$countTable,delim="\t",col_names = T,skip = 0)
print("input countTable:")
[1] "input countTable:"
print(dim(t))
[1] 22246     5
t[1:3,1:3]
# A tibble: 3 × 3
  sites  `MORC3-V5_DMSO_M3AID_REP1_featureCount` `MORC3-V5_DMSO_M3AID_REP2_featureCount`
  <chr>                                    <dbl>                                   <dbl>
1 site_1                                      44                                      26
2 site_2                                      72                                      47
3 site_3                                      26                                      13
t[t$sites=="site_1468",]
# A tibble: 1 × 5
  sites     `MORC3-V5_DMSO_M3AID_REP1_featureCount` `MORC3-V5_DMSO_M3AID_REP2_featureCount` sgAAVS1_DMSO_R1_featureCount sgAAVS1_DMSO_R2_featureCount
  <chr>                                       <dbl>                                   <dbl>                        <dbl>                        <dbl>
1 site_1468                                      29                                      14                           10                           22
#
# sample meta info
#
m <- read_delim(params$sampleTable,delim="\t",col_names = T)
m
# A tibble: 4 × 4
  sample_id sample_name                   count_column                          condition
      <dbl> <chr>                         <chr>                                 <chr>    
1         1 MORC3_V5_DMSO_M3AID_REP1_exp2 MORC3-V5_DMSO_M3AID_REP1_featureCount exp2     
2         2 MORC3_V5_DMSO_M3AID_REP2_exp2 MORC3-V5_DMSO_M3AID_REP2_featureCount exp2     
3         3 sgAAVS1_DMSO_R1_exp1          sgAAVS1_DMSO_R1_featureCount          exp1     
4         4 sgAAVS1_DMSO_R2_exp1          sgAAVS1_DMSO_R2_featureCount          exp1     
datatable(m)
#
# annotation data
#
countTableAnnot <- read_delim(params$countTableAnnot,delim="\t",col_names = T)
colnames(countTableAnnot)[1] <- "SiteID"
countTableAnnot[1:3,]
# A tibble: 3 × 6
  SiteID Chr    Start    End Strand Length
  <chr>  <chr>  <dbl>  <dbl> <chr>   <dbl>
1 site_1 chr1  826913 827797 +         885
2 site_2 chr1  904446 905303 +         858
3 site_3 chr1  906294 907109 +         816
countTableAnnot %>% filter(SiteID=="site_1468")
# A tibble: 1 × 6
  SiteID    Chr       Start       End Strand Length
  <chr>     <chr>     <dbl>     <dbl> <chr>   <dbl>
1 site_1468 chr1  146472327 146472590 +         264
# counts only
tc <- t[,-1]
print(dim(tc))
[1] 22246     4
# setup DESeqDataSet
tcm <- as.matrix(tc)
tcm[1:3,1:3]
     MORC3-V5_DMSO_M3AID_REP1_featureCount MORC3-V5_DMSO_M3AID_REP2_featureCount sgAAVS1_DMSO_R1_featureCount
[1,]                                    44                                    26                           60
[2,]                                    72                                    47                           88
[3,]                                    26                                    13                           76
counts <- matrix(tcm,
                 ncol=dim(tcm)[2],
                 dimnames=list(t$sites,colnames(tcm)))
print("Final deseq2 count table:")
[1] "Final deseq2 count table:"
print(dim(counts))
[1] 22246     4
counts[1:3,]
       MORC3-V5_DMSO_M3AID_REP1_featureCount MORC3-V5_DMSO_M3AID_REP2_featureCount sgAAVS1_DMSO_R1_featureCount sgAAVS1_DMSO_R2_featureCount
site_1                                    44                                    26                           60                           69
site_2                                    72                                    47                           88                           88
site_3                                    26                                    13                           76                           39

2 DESeq2

DDS <- DESeqDataSetFromMatrix(countData = counts, 
                             colData = m,
                             design = ~condition)
head(colData(DDS))
DataFrame with 4 rows and 4 columns
                                      sample_id            sample_name           count_column condition
                                      <numeric>            <character>            <character>  <factor>
MORC3-V5_DMSO_M3AID_REP1_featureCount         1 MORC3_V5_DMSO_M3AID_.. MORC3-V5_DMSO_M3AID_..      exp2
MORC3-V5_DMSO_M3AID_REP2_featureCount         2 MORC3_V5_DMSO_M3AID_.. MORC3-V5_DMSO_M3AID_..      exp2
sgAAVS1_DMSO_R1_featureCount                  3   sgAAVS1_DMSO_R1_exp1 sgAAVS1_DMSO_R1_feat..      exp1
sgAAVS1_DMSO_R2_featureCount                  4   sgAAVS1_DMSO_R2_exp1 sgAAVS1_DMSO_R2_feat..      exp1
rowData(DDS) <- countTableAnnot
head(rowData(DDS))
DataFrame with 6 rows and 6 columns
            SiteID         Chr     Start       End      Strand    Length
       <character> <character> <numeric> <numeric> <character> <numeric>
site_1      site_1        chr1    826913    827797           +       885
site_2      site_2        chr1    904446    905303           +       858
site_3      site_3        chr1    906294    907109           +       816
site_4      site_4        chr1    921152    921458           +       307
site_5      site_5        chr1    924657    924957           +       301
site_6      site_6        chr1    940863    941872           +      1010
print(dim(DDS))
[1] 22246     4
rowData(DDS)[rowData(DDS)$SiteID=="site_1468",]
DataFrame with 1 row and 6 columns
               SiteID         Chr     Start       End      Strand    Length
          <character> <character> <numeric> <numeric> <character> <numeric>
site_1468   site_1468        chr1 146472327 146472590           +       264

2.1 Relative fold changes and PCA

The Relative Log Expression (RLE) plot is a useful diagnostic plot to visualize the differences between the distributions of read counts across samples.

It shows the boxplots of the log-ratios of the site-level read counts of each sample to those of a reference sample (defined as the median across the samples). Ideally, the distributions should be centered around the zero line and as tight as possible. Clear deviations indicate the need for normalization and/or the presence of outlying samples.

The PCA plot represents a Principal Component Analysis (PCA) plot of the raw counts in object .

#
# count check
#
par(las=2,mar=c(10,2,2,2))
plotRLE(counts(DDS))

plotPCA(counts(DDS))

#
# count distribution
#

# all counts of all samples together
raw_count_cuttoff <- 0

as_tibble(counts(DDS)) %>% 
  pivot_longer(cols=everything(),names_to = "samples", values_to = "counts") %>%
  ggplot(aes(x=log2(counts+1))) +
  geom_histogram(bins=100) +
  geom_vline(xintercept = log2(raw_count_cuttoff),col="red") +
  theme_light() +
  ggtitle(paste("raw count distributions (all sites), cutoff=",raw_count_cuttoff))

raw_count_cuttoff <- 1

as_tibble(counts(DDS)) %>% 
  pivot_longer(cols=everything(),names_to = "samples", values_to = "counts") %>%
  ggplot(aes(x=log2(counts+1))) +
  geom_histogram(bins=100) +
  geom_vline(xintercept = log2(raw_count_cuttoff),col="red") +
  theme_light() +
  ggtitle(paste("raw count distributions (all sites), cutoff=",raw_count_cuttoff))

# by sample
as_tibble(counts(DDS)) %>% 
  pivot_longer(cols=everything(),names_to = "samples", values_to = "counts") %>%
  ggplot(aes(x=log2(counts+1))) +
  geom_histogram(bins=100) +
  geom_vline(xintercept = log2(raw_count_cuttoff),col="red") +
  facet_wrap(~samples) +
  theme_light() +
  ggtitle(paste("raw count distributions (all sites), cutoff=",raw_count_cuttoff))

# rowsums
rs <- tibble(raw_count_sums_per_peak=rowSums(counts(DDS)))
rs %>% ggplot(aes(x=log2(raw_count_sums_per_peak+1))) + 
  geom_histogram(bins=100) +
  theme_light() +
  ggtitle("Raw count rowsums (over all replicates, per site)")

#
# filter a bit for very low numbers
#
print("low counts filtered out:")
[1] "low counts filtered out:"
smallestGroupSize <- 2

keep <- rowSums(counts(DDS) > (2*raw_count_cuttoff)) >= smallestGroupSize
print(table(keep))
keep
FALSE  TRUE 
    2 22244 
DDS <- DDS[keep,]
dim(DDS)
[1] 22244     4
plotRLE(counts(DDS))

plotPCA(counts(DDS))

sections=c(10,40,800,2000,5000)
plot_counts(counts(DDS)+1,sections)
[1] "plotting bargraph ..."
[1] "plotting overall count sum ditribution ..."
[1] "plotting overall count sum ditribution, by time and replicate ..."
[1] "final plot ..."

2.2 RLOG transform (VST)

#
# get rlogs
#

if(file.exists(paste0(params$odir,"/vst.RData"))) {
  print(paste0("loading ",params$odir,"/vst.RData"))
  load(paste0(params$odir,"/vst.RData"))
} else {
  print("calculating VST ...")
  vst <- rlog(DDS, blind=T)
  print(dim(DDS))
  print(dim(vst))
  save(vst,file=paste0(params$odir,"/vst.RData"))
}
[1] "calculating VST ..."
[1] 22244     4
[1] 22244     4

2.2.1 PCA

We investigate the relationship between mean and variance, and its effect on the PCA.

time.var.png

mPCA(m=assay(vst),
  colgroups=m$condition,
  collegend="Condition",
  pchgroups=as.character(m$sample_id),
  pchlegend="Replicate",
  #annot=rowData(vst)[,1],
  title="vst.blind",
  ofile=paste0(params$odir,"/",params$condition),
  varCutoff = 0.9,
  show_biplot=F)
[1] "mPCA ..."
[1] "make matrix ..."
[1] 22244     4
[1] "pchgroups set to:"
[1] "1" "2" "3" "4"
[1] "pchgroups set to:"
1 2 3 4 
1 2 3 4 
[1] "calculate variance ..."
[1] "plotting variances ..."
[1] "violin ..."
[1] "density ..."
[1] "plotting variances versus means ..."
[1] "density ..."
[1] "scatter ..."
[1] "patchwork ..."
[1] "selected features (variances only):"
[1] 2225
[1] "top selected features (variances only):"
NULL
[1] "taking default labels ...."
[1] "labels:"
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount" "MORC3-V5_DMSO_M3AID_REP2_featureCount" "sgAAVS1_DMSO_R1_featureCount"          "sgAAVS1_DMSO_R2_featureCount"         
[1] "pca ..."
[1] "plotIndiv ..."
[1] "colgroups_factor:"
[1] "exp2" "exp2" "exp1" "exp1"
[1] "pchgroup:"
1 2 3 4 
1 2 3 4 
[1] "pchlevels:"
[1] "1" "2" "3" "4"

[1] "selected features (variance and means):"
[1] "0.0703069724604185 5.35432301347482"
[1] 1214
[1] "top selected features (variance and means):"
NULL

[1] "selected features (no filtering, only variance == 0 filtered out):"
[1] 22244
[1] "top selected features (no filtering, only variance == 0 filtered out):"
NULL
[1] "selected features (variance>0):"
[1] 22244

2.3 Pairwise scatters RAW & VST

#
# plot for each celltype
#

# subsample, for quicker drawing
if(params$downsample_pairs>0) {
  si <- sample(c(1:dim(DDS)[1]),params$downsample_pairs)
  print(head(si))
  DDS.sub <- DDS[si,]
} else {
  DDS.sub <- DDS
  si <- c(1:dim(DDS)[1])
}
[1] 18753 21657  9290  1252 15506  8826
# plot raw counts
p <- ggpairs(data.frame(log2(counts(DDS.sub)+1)),
      lower=list(continuous=my_fn)) +
  ggtitle(paste0("Raw counts subsample: ",params$downsample_pairs)) +
  theme_bw()
#p_(p)
print(p)

ggsave(p,filename=paste0(params$odir,"/raw.ss_",params$downsample_pairs,".jpeg")) 
  
# plot vst counts
vst.sub <- vst[si,]
p <- ggpairs(data.frame(assay(vst.sub)),
      lower=list(continuous=my_fn)) +
ggtitle(paste0("RLOG transformed counts, subsample: ",params$downsample_pairs)) +
  theme_bw()
#p_(p)
print(p)

ggsave(p,filename=paste0(params$odir,"/vst.ss_",params$downsample_pairs,".jpeg")) 

2.4 Pairwise scatters TPM

#
# plot for each celltype
#

# subsample, for quicker drawing
if(params$downsample_pairs>0) {
  si <- sample(c(1:dim(DDS)[1]),params$downsample_pairs)
  print(head(si))
  DDS.sub <- DDS[si,]
} else {
  DDS.sub <- DDS
  si <- c(1:dim(DDS)[1])
}
[1]  7580  6620  4288 16478  5322 14895
# plot TPMs
rd <- rowData(DDS.sub)
sitelengths <- rd$End-rd$Start
tpm.counts <- convertCounts(countsMatrix = counts(DDS.sub),unit="TPM",geneLength=sitelengths)

p <- ggpairs(data.frame(log2(tpm.counts+1)),
      lower=list(continuous=my_fn)) +
  ggtitle(paste0("TPMs, subsample: ",params$downsample_pairs)) +
  theme_bw()
#p_(p)
print(p)

ggsave(p,filename=paste0(params$odir,"/tpm.ss_",params$downsample_pairs,".jpeg")) 

2.5 Dispersion and normalized counts

plotDispEsts plots the per-gene dispersion estimates together with the fitted mean-dispersion relationship.

#
# ---- DESeq2 model ----
#
#if(file.exists(paste0(params$odir,"/deobj.",params$celltype,",RData"))) {
#  print(paste0("loading ",params$odir,"/deobj.",params$celltype,".RData"))
#  load(paste0(params$odir,"/deobj.",params$celltype,".RData"))
#} else {
  print("calculating DESeq2 object ...")
[1] "calculating DESeq2 object ..."
  deobj <- DESeq(DDS, test='Wald')
  save(deobj,file=paste0(params$odir,"/deobj.RData"))
#}

print("DESeq2 sizeFactors:")
[1] "DESeq2 sizeFactors:"
p <- tibble("names"=names(sizeFactors(deobj)),
     "sf"=sizeFactors(deobj)) %>%
 ggplot(aes(x=names,y=sf)) +
 geom_col() +
 theme_light() +
 theme(axis.text.x=element_text(angle=90)) +
 ggtitle("Size factors")
print(p)

ggsave(p,filename=paste0(params$odir,"/sizefactors.jpeg")) 
# 
print("Dispersion estimates:")
[1] "Dispersion estimates:"
plotDispEsts(deobj)

res <- results(deobj)

2.6 Relative fold change and PCA - normalized

par(las=2,mar=c(10,2,2,2))
plotRLE(counts(deobj,normalized=T))

plotPCA(counts(deobj,normalized=T))

2.7 Pairwise scatters NORM

#
# ---- DESeq2 plots ----
#

# ---- norm counts ----
cts.norm <- counts(deobj,normalized=T)

# same si as above
if(params$downsample_pairs>0) {
  cts.norm.sub <- cts.norm[si,]
} else {
  cts.norm.sub <- cts.norm
}
# summary(is.na(cts.norm)) # there should be no NA
p <- ggpairs(log2(cts.norm.sub+1),
      lower=list(continuous=my_fn)) +
  ggtitle(paste0("Normalized counts, subsample=",params$downsample_pairs)) +
  theme_bw()
#p_(p)
print(p)

ggsave(p,filename=paste0(params$odir,"/norm.ss_",params$downsample_pairs,".jpeg")) 

2.8 MA and volcano plots

if(params$downsample_ma>0) {
  si <- sample(c(1:dim(deobj)[1]),params$downsample_ma)
  deobj.sub <- deobj[si,]
  vst.sub <- vst[si,]
} else {
  si <- c(1:dim(deobj)[1])
  deobj.sub <- deobj
  vst.sub <- vst
}

res <- results(deobj.sub,contrast = c("condition",params$experiment_condition,params$control_condition))
res.df <- results(deobj.sub,contrast = c("condition",params$experiment_condition,params$control_condition),tidy=T)

dim(res.df)
[1] 22244     7
res.df[1:3,1:3]
     row baseMean log2FoldChange
1 site_1 46.75372     0.04882633
2 site_2 71.95839     0.37476977
3 site_3 33.45225    -0.62945101
res.df[1:3,]
     row baseMean log2FoldChange     lfcSE        stat    pvalue      padj
1 site_1 46.75372     0.04882633 0.5478855  0.08911775 0.9289883 0.9994137
2 site_2 71.95839     0.37476977 0.4462411  0.83983702 0.4009998 0.9994137
3 site_3 33.45225    -0.62945101 0.6765031 -0.93044807 0.3521391        NA
#data.frame(rowData(deobj.sub))
row.df <- data.frame(rowData(deobj.sub)[1:19])
row.df <- cbind(row=rownames(row.df),row.df)
row.df[1:3,1:3]
          row SiteID  Chr
site_1 site_1 site_1 chr1
site_2 site_2 site_2 chr1
site_3 site_3 site_3 chr1
res.df.a <- merge(x=row.df,y=res.df,by="row",sort=T,all=T)
res.df.a[1:3,]
       row   SiteID  Chr   Start     End Strand Length baseMean.x   baseVar allZero dispGeneEst dispGeneIter   dispFit dispersion dispIter dispOutlier   dispMAP Intercept condition_exp2_vs_exp1 SE_Intercept baseMean.y
1   site_1   site_1 chr1  826913  827797      +    885   46.75372  39.33170   FALSE  0.00000001            2 0.1462236  0.1218074        6       FALSE 0.1218074  5.524677             0.04882633    0.3779915   46.75372
2  site_10  site_10 chr1  999572 1000802      +   1231   50.08206 454.39487   FALSE  0.26045883            3 0.1373907  0.1567706        9       FALSE 0.1567706  5.755564            -0.21399560    0.4205713   50.08206
3 site_100 site_100 chr1 2953145 2953475      +    331   10.95914  80.33834   FALSE  1.81124086            7 0.5803305  0.6815258        9       FALSE 0.6815258  3.980974            -1.32309113    0.8696165   10.95914
  log2FoldChange     lfcSE        stat    pvalue      padj
1     0.04882633 0.5478855  0.08911775 0.9289883 0.9994137
2    -0.21399560 0.6093170 -0.35120570 0.7254340 0.9994137
3    -1.32309113 1.3021971 -1.01604523 0.3096078        NA
dim(res.df.a)
[1] 22244    26
dim(vst.sub)
[1] 22244     4
dim(counts(deobj.sub))
[1] 22244     4
desc <- "test"
write_extended_results_table(deobj=deobj.sub,
                              vst=vst.sub,
                              res=res.df.a,
                              ofile=paste0(params$odir,"/",desc,".",params$downsample_ma,".txt"),
                             genelength = res.df.a$End-res.df.a$Start)
[1] 22244
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
[1] 22244     4
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
[1] 22244     4
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
[1] 22244     4
[1] "MORC3-V5_DMSO_M3AID_REP1_featureCount"
[1] "MORC3-V5_DMSO_M3AID_REP2_featureCount"
[1] "sgAAVS1_DMSO_R1_featureCount"
[1] "sgAAVS1_DMSO_R2_featureCount"
#
# MA plots
#
res
log2 fold change (MLE): condition exp2 vs exp1 
Wald test p-value: condition exp2 vs exp1 
DataFrame with 22244 rows and 6 columns
            baseMean log2FoldChange     lfcSE        stat    pvalue      padj
           <numeric>      <numeric> <numeric>   <numeric> <numeric> <numeric>
site_1       46.7537     0.04882633  0.547886  0.08911775  0.928988  0.999414
site_2       71.9584     0.37476977  0.446241  0.83983702  0.401000  0.999414
site_3       33.4523    -0.62945101  0.676503 -0.93044807  0.352139        NA
site_4       21.0182    -0.00537712  0.801507 -0.00670876  0.994647        NA
site_5       16.1062    -0.67071868  0.949427 -0.70644553  0.479911        NA
...              ...            ...       ...         ...       ...       ...
site_22242   19.6523      -0.331930  0.846293   -0.392217  0.694898        NA
site_22243   17.4983      -1.117983  0.888502   -1.258279  0.208291        NA
site_22244   49.2172      -0.271975  0.615693   -0.441737  0.658679  0.999414
site_22245   19.8220      -0.191921  0.826885   -0.232101  0.816460        NA
site_22246   30.5679       0.191471  0.702078    0.272720  0.785069        NA
p <- my_ma(padj_co=0.01,deobj=res)

print(p)
NULL
p1 <- my_ma(padj_co=0.01,deobj=NULL,df=data.frame(res),baseMean_co=0, color=T,
            filter_df_na = T, name=desc,
            outfile=paste0(params$odir,"/ma.",params$downsample_ma,".jpeg"))
print(p1)

p2 <- my_volcano(df=data.frame(res),baseMean_co = 0,filter_df_na = T, 
                 color=T, name="test", 
                 outfile=paste0(params$odir,"/vc.",params$downsample_ma,".jpeg"))
print(p2)

# top results table with links
session="http://localhost:60151/goto?"
pad=10000
html.df <- coords_to_igv_links(df = res.df.a %>% arrange(pvalue) %>% slice_head(n=100),
                                chr="Chr",start="Start",end="End",session=session,pad=pad) # params$session)
html.table <- html.df %>% 
  addHtmlTableStyle(align = "r") %>%
  htmlTable(caption=paste("test","top 100 entries, sorted by pvalue, pad=",pad),
            useViewer=FALSE,rnames=FALSE)
write(html.table,file = paste0(params$odir,"/test.",params$downsample_ma,".top100.html"))

3 Session info

sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS;  LAPACK version 3.11.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Vienna
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DGEobj.utils_1.0.6          htmlTable_2.4.2             DT_0.31                     data.table_1.15.4           ggpointdensity_0.1.0        patchwork_1.2.0             viridis_0.6.5              
 [8] viridisLite_0.4.2           RColorBrewer_1.1-3          mixOmics_6.26.0             lattice_0.22-6              MASS_7.3-60.0.1             EDASeq_2.36.0               apeglm_1.24.0              
[15] ShortRead_1.60.0            GenomicAlignments_1.38.2    Rsamtools_2.18.0            Biostrings_2.70.3           XVector_0.42.0              BiocParallel_1.36.0         ape_5.8                    
[22] ReportingTools_2.42.3       knitr_1.45                  vsn_3.70.0                  edgeR_4.0.16                limma_3.58.1                DESeq2_1.42.1               SummarizedExperiment_1.32.0
[29] Biobase_2.62.0              MatrixGenerics_1.14.0       matrixStats_1.3.0           GenomicRanges_1.54.1        GenomeInfoDb_1.38.8         IRanges_2.36.0              S4Vectors_0.40.2           
[36] BiocGenerics_0.48.1         ggrepel_0.9.5               GGally_2.2.1                lubridate_1.9.3             forcats_1.0.0               stringr_1.5.0               dplyr_1.1.4                
[43] purrr_1.0.2                 readr_2.1.5                 tidyr_1.3.1                 tibble_3.2.1                ggplot2_3.5.1               tidyverse_2.0.0            

loaded via a namespace (and not attached):
  [1] ProtGenerics_1.34.0      bitops_1.0-7             httr_1.4.7               Rgraphviz_2.46.0         numDeriv_2016.8-1.1      tools_4.3.2              backports_1.4.1          utf8_1.2.4              
  [9] R6_2.5.1                 mgcv_1.9-1               lazyeval_0.2.2           withr_3.0.0              prettyunits_1.2.0        gridExtra_2.3            preprocessCore_1.64.0    textshaping_0.3.7       
 [17] cli_3.6.2                labeling_0.4.3           ggbio_1.50.0             sass_0.4.7               mvtnorm_1.2-4            genefilter_1.84.0        systemfonts_1.0.5        foreign_0.8-86          
 [25] R.utils_2.12.3           AnnotationForge_1.44.0   dichromat_2.0-0.1        BSgenome_1.70.2          bbmle_1.0.25.1           rstudioapi_0.15.0        RSQLite_2.3.6            generics_0.1.3          
 [33] GOstats_2.68.0           BiocIO_1.12.0            hwriter_1.3.2.1          crosstalk_1.2.1          vroom_1.6.5              GO.db_3.18.0             Matrix_1.6-5             interp_1.1-6            
 [41] fansi_1.0.5              abind_1.4-5              R.methodsS3_1.8.2        lifecycle_1.0.4          yaml_2.3.7               SparseArray_1.2.4        BiocFileCache_2.10.2     grid_4.3.2              
 [49] blob_1.2.4               crayon_1.5.2             bdsmatrix_1.3-7          GenomicFeatures_1.54.4   annotate_1.80.0          KEGGREST_1.42.0          pillar_1.9.0             rjson_0.2.21            
 [57] corpcor_1.6.10           codetools_0.2-20         glue_1.7.0               vctrs_0.6.5              png_0.1-8                gtable_0.3.5             assertthat_0.2.1         emdbook_1.3.13          
 [65] cachem_1.0.8             aroma.light_3.32.0       xfun_0.41                S4Arrays_1.2.1           coda_0.19-4.1            survival_3.6-4           statmod_1.5.0            ellipsis_0.3.2          
 [73] nlme_3.1-164             Category_2.68.0          bit64_4.0.5              progress_1.2.3           filelock_1.0.3           bslib_0.5.1              affyio_1.72.0            rpart_4.1.23            
 [81] colorspace_2.1-0         DBI_1.2.2                Hmisc_5.1-2              nnet_7.3-19              tidyselect_1.2.1         bit_4.0.5                compiler_4.3.2           curl_5.1.0              
 [89] graph_1.80.0             xml2_1.3.5               DelayedArray_0.28.0      rtracklayer_1.62.0       checkmate_2.3.1          scales_1.3.0             hexbin_1.28.3            affy_1.80.0             
 [97] RBGL_1.78.0              rappdirs_0.3.3           digest_0.6.33            rmarkdown_2.25           htmltools_0.5.7          pkgconfig_2.0.3          jpeg_0.1-10              base64enc_0.1-3         
[105] highr_0.10               dbplyr_2.5.0             fastmap_1.1.1            ensembldb_2.26.0         rlang_1.1.3              htmlwidgets_1.6.2        PFAM.db_3.18.0           farver_2.1.1            
[113] jquerylib_0.1.4          jsonlite_1.8.7           R.oo_1.26.0              VariantAnnotation_1.48.1 RCurl_1.98-1.14          magrittr_2.0.3           Formula_1.2-5            GenomeInfoDbData_1.2.11 
[121] munsell_0.5.1            Rcpp_1.0.11              stringi_1.7.12           zlibbioc_1.48.2          plyr_1.8.9               ggstats_0.6.0            parallel_4.3.2           deldir_2.0-4            
[129] splines_4.3.2            hms_1.1.3                locfit_1.5-9.9           igraph_2.0.3             reshape2_1.4.4           biomaRt_2.58.2           DGEobj_1.1.2             XML_3.99-0.16.1         
[137] evaluate_0.23            latticeExtra_0.6-30      biovizBase_1.50.0        BiocManager_1.30.22      tzdb_0.4.0               xtable_1.8-4             restfulr_0.0.15          AnnotationFilter_1.26.0 
[145] RSpectra_0.16-1          ragg_1.2.6               rARPACK_0.11-0           OrganismDbi_1.44.0       memoise_2.0.1            AnnotationDbi_1.64.1     ellipse_0.5.0            cluster_2.1.6           
[153] timechange_0.3.0         GSEABase_1.64.0